set more off


* the data used in this paper may not be publicly released and aspects of the data cleaning files may violate terms of our data use agreeements.
* our code takes a final analysis dataset as given, and provides code of all actual analyses in the paper.
* this file reflects analysis for period 2 (2006-2012), bottled water, batteries, flashlights.

**************************************************************************************************
* LEVEL BY WEEK ANALYSIS CONSTRUCTION 2006-2012
**************************************************************************************************

	**************************************************************************************************
	* THIS IS WHERE LEVEL IS CHOSEN
	**************************************************************************************************
	
* level is county (0) or store (1)
local level=1

use period2_county_treatments, clear 
tsset, clear

if `level'==1 {
	merge 1:m countyid statadate using period2_store_sales	 
}
else {
	merge 1:1 countyid statadate using period2_county_sales
}

keep if _merge==3
	* merge=2 in both datasets is weeks without full info ... jan 2006, dec 2012.
	* merge=2 in store-level data is virginia cities without threat information
	* merge=1 is counties without stores that fit criteria
drop _merge

destring countyid, gen(countyid_num)

sort statadate
egen week=group(statadate)

if `level' == 1 {
	tsset store_code_uc week	 
} 
else {
	tsset countyid_num week
}


******************************************************************************************************
* EVENT STUDY
*****************************************************************************************************

* note that these are point-wise confidence intervals

if `level'==1 { 

foreach z in "1487" {

preserve

sum expend_`z'
gen mean=r(mean)
replace expend_`z'=expend_`z'-mean
replace expend_`z'=expend_`z'/1000

local win = 10
forvalues x=0/`win' {
	gen expend_`z'_l`x'=l`x'.expend_`z'
}
forvalues x=1/`win' {
	gen expend_`z'_f`x'=f`x'.expend_`z'
}

*keep if hlandfall_false==1
keep if hlandfall==1

keep expend_`z'_*

forvalues x=0/`win' {
	gen expend2_`z'_l`x'=expend_`z'_l`x'
}
forvalues x=1/`win' {
	gen expend2_`z'_f`x'=expend_`z'_f`x'
}

collapse (sem) expend2_`z'_* (mean) expend_`z'_*

duplicates drop
forvalues x=1/`win' {
	local y=`x'+50
	rename expend_`z'_f`x' expend_`z'_`y'
	rename expend2_`z'_f`x' expend2_`z'_`y'
}
forvalues x=0/`win' {
	local y=50-`x'
	rename expend_`z'_l`x' expend_`z'_`y'
	rename expend2_`z'_l`x' expend2_`z'_`y'
}

gen id=1
reshape long expend_`z'_ expend2_`z'_, i(id) j(time)
drop id
replace time=time-50

gen ci=expend2_`z'*1.96
gen ul=expend_`z'+ci
gen ll=expend_`z'-ci

twoway (rarea ul ll time, color(gs14)) (connected expend_`z'_ time, msize(small)) (pcarrowi 0 0 10 0, msize(zero)), scheme(s2mono) plotr(m(zero)) ytitle("") ylabel(0(3)9) xlabel(-10(2)10) xtitle("Weeks Before and After Storm Landfall") legend(off) title("Water Sales - 2006-2012") saving(Fig2B, replace) 

restore

}


foreach z in "7870" {

preserve

sum expend_`z'
gen mean=r(mean)
replace expend_`z'=expend_`z'-mean
replace expend_`z'=expend_`z'/1000

local win = 10
forvalues x=0/`win' {
	gen expend_`z'_l`x'=l`x'.expend_`z'
}
forvalues x=1/`win' {
	gen expend_`z'_f`x'=f`x'.expend_`z'
}

*keep if hlandfall_false==1
keep if hlandfall==1

keep expend_`z'_*

forvalues x=0/`win' {
	gen expend2_`z'_l`x'=expend_`z'_l`x'
}
forvalues x=1/`win' {
	gen expend2_`z'_f`x'=expend_`z'_f`x'
}

collapse (sem) expend2_`z'_* (mean) expend_`z'_*

duplicates drop
forvalues x=1/`win' {
	local y=`x'+50
	rename expend_`z'_f`x' expend_`z'_`y'
	rename expend2_`z'_f`x' expend2_`z'_`y'
}
forvalues x=0/`win' {
	local y=50-`x'
	rename expend_`z'_l`x' expend_`z'_`y'
	rename expend2_`z'_l`x' expend2_`z'_`y'
}

gen id=1
reshape long expend_`z'_ expend2_`z'_, i(id) j(time)
drop id
replace time=time-50

gen ci=expend2_`z'*1.96
gen ul=expend_`z'+ci
gen ll=expend_`z'-ci

twoway (rarea ul ll time, color(gs14)) (connected expend_`z'_ time, msize(small)) (pcarrowi 0 0 3 0, msize(zero)), scheme(s2mono) plotr(m(zero)) ytitle("") ylabel(0(1)3) xlabel(-10(2)10) xtitle("Weeks Before and After Storm Landfall") legend(off) title("Battery Sales - 2006-2012") saving(Fig2C, replace)

restore

}


foreach z in "7871" {

preserve

sum expend_`z'
gen mean=r(mean)
replace expend_`z'=expend_`z'-mean
replace expend_`z'=expend_`z'/1000

local win = 10
forvalues x=0/`win' {
	gen expend_`z'_l`x'=l`x'.expend_`z'
}
forvalues x=1/`win' {
	gen expend_`z'_f`x'=f`x'.expend_`z'
}

*keep if hlandfall_false==1
keep if hlandfall==1

keep expend_`z'_*

forvalues x=0/`win' {
	gen expend2_`z'_l`x'=expend_`z'_l`x'
}
forvalues x=1/`win' {
	gen expend2_`z'_f`x'=expend_`z'_f`x'
}

collapse (sem) expend2_`z'_* (mean) expend_`z'_*

duplicates drop
forvalues x=1/`win' {
	local y=`x'+50
	rename expend_`z'_f`x' expend_`z'_`y'
	rename expend2_`z'_f`x' expend2_`z'_`y'
}
forvalues x=0/`win' {
	local y=50-`x'
	rename expend_`z'_l`x' expend_`z'_`y'
	rename expend2_`z'_l`x' expend2_`z'_`y'
}

gen id=1
reshape long expend_`z'_ expend2_`z'_, i(id) j(time)
drop id
replace time=time-50

gen ci=expend2_`z'*1.96
gen ul=expend_`z'+ci
gen ll=expend_`z'-ci

twoway (rarea ul ll time, color(gs14)) (connected expend_`z'_ time, msize(small)) (pcarrowi 0 0 0.3 0, msize(zero)), scheme(s2mono) plotr(m(zero)) ytitle("") ylabel(0(.1).3) xlabel(-10(2)10) xtitle("Weeks Before and After Storm Landfall") legend(off) title("Flashlight Sales - 2006-2012") saving(Fig2D, replace)

restore

}

}



****************************************************************************************************
* ADDITIONAL PROCESSING FOR REGRESSION ANALYSIS
****************************************************************************************************


* merge in weather at county level
* some winter weeks are missing weather data. we do not use winter weeks, so do not worry about this 
* merge==2 is  weeks with no sales
if `level' == 1 {
	merge m:1 countyid statadate using period2_county_weather	 
} 
else {
	merge 1:1 countyid statadate using period2_county_weather
}
drop if _merge==2
drop _merge

* create year dummies
tab year, gen(years)

* generate logged depvars - NOTE ADD $1 and log
foreach code of numlist 1487 7460 7870 7871 8420 8301 6008 6012 5502 {
	gen lexpend_`code'=log(expend_`code'+1)
}

foreach stub in "w16906" "w16912" "w128001" "b2ct" "b4ct" "b8ct" "b16ct" "f1ct" "f2ct" "f3ct" "f4ct" {
	gen lexpend_`stub'=log(expend_`stub'+1)
}

* create month and day varaibles
capture drop month
capture drop day
gen month=month(statadate)
gen day=day(statadate)

* define first week of hurricane season 
if `level' == 0 {

gen first=((year==2006&month==5&day==28)|(year==2007&month==5&day==27)|(year==2008&month==6&day==1)|(year==2009&month==5&day==31)|(year==2010&month==5&day==30)|(year==2011&month==5&day==29)|(year==2012&month==5&day==27))

* define week of first hurricane of that year's hurricane season. 
* set to missing if that county is impacted by a hurricane sometime during that year. 
	bys year: egen temp=min(statadate) if hlandfall==1
	bys year: egen temp2 =min(temp)
	format temp2 %td
	gen firsth=(temp2==statadate)
	drop temp temp2
	bys countyid year: egen temp1=max(htreat_before)
	bys countyid year: egen temp2=max(htreat_after)
	bys countyid year: egen temp3=max(hlandfall)
	egen ever=rowmax(temp1 temp2 temp3)
	replace firsth=. if ever==1
	drop temp1 temp2 temp3

}


* merge in demos at county level
merge m:1 countyid using county_demos
keep if _m==3
drop _merge

* merge in risk factors at county level 
merge m:1 countyid using county_coastal
drop if _merge==2
gen coastal_county=(coastal=="C")
drop coastal _merge

if `level' == 1 {
	merge m:1 countyid statadate using period2_county_fema	 
} 
else {
	merge 1:1 countyid statadate using period2_county_fema
}
drop if _merge==2
gen fema_county_week=(fema==1)
drop fema _merge

rename pctbachelorsormore pctbach


foreach var in pctblack1 medianhhinc pctbach {
	quietly sum `var', detail
	scalar m`var'=r(p50)
	gen hi`var'=(`var'>m`var') if `var'!=.
	gen b_inter_hi`var'=hi`var'*htreat_before
	gen d_inter_hi`var'=hi`var'*hlandfall
	gen a_inter_hi`var'=hi`var'*htreat_after
}

if `level' == 0 {

	foreach var in pctblack1 medianhhinc pctbach {
		gen first_inter_hi`var'=hi`var'*first
		gen firsth_inter_hi`var'=hi`var'*firsth
	}
}


foreach rvar in coastal_county fema_county_week {
	gen b_inter_`rvar'=`rvar'*htreat_before
	gen d_inter_`rvar'=`rvar'*hlandfall
	gen a_inter_`rvar'=`rvar'*htreat_after
	gen b_inter_`rvar'_evac=`rvar'*htreat_before_evac
	gen d_inter_`rvar'_evac=`rvar'*hlandfall_evac
	gen a_inter_`rvar'_evac=`rvar'*htreat_after_evac
	gen b_inter_`rvar'_nonevac=`rvar'*htreat_before_nonevac
	gen d_inter_`rvar'_nonevac=`rvar'*hlandfall_nonevac
	gen a_inter_`rvar'_nonevac=`rvar'*htreat_after_nonevac
}

if `level' == 0 {

	foreach rvar in coastal_county fema_county_week {
		gen first_inter_`rvar'=`rvar'*first
		gen firsth_inter_`rvar'=`rvar'*firsth
	}
}

* assign fixed effects level
if `level'==1 {
	local FE "store_code_uc"	 
}
else {
	local FE "countyid_num"
}

* generate impact leads
if `level' == 1 {
	tsset store_code_uc week	 
} 
else {
	tsset countyid_num week
}

gen hlandfall_f1=f1.hlandfall
gen hlandfall_f2=f2.hlandfall


***********************************************************************************************************************************
* REGRESSIONS
***********************************************************************************************************************************

* early season stockpiling 

if `level' == 0 {

preserve

	* this analysis must be done before full restrictions to hurricane season only.
	* here, we restrict analysis to april - november to allow some pre-season data but to otherwise keep analysis comparable.

keep if inrange(month,5,11)
tab month, gen(months)

foreach code in "1487" "7870" "7871" {
	xtreg lexpend_`code' first prcp tmean tmax tmin months1-months6 years1-years6,  fe i(`FE') vce(cluster countyid_num)
	xtreg lexpend_`code' first first_inter_hipctblack1 prcp tmean tmax tmin months1-months6 years1-years6,  fe i(`FE') vce(cluster countyid_num)
	xtreg lexpend_`code' first first_inter_himedianhhinc prcp tmean tmax tmin months1-months6 years1-years6,  fe i(`FE') vce(cluster countyid_num)
	xtreg lexpend_`code' first first_inter_hipctbach prcp tmean tmax tmin months1-months6 years1-years6,  fe i(`FE') vce(cluster countyid_num)
	xtreg lexpend_`code' first first_inter_coastal_county prcp tmean tmax tmin months1-months6 years1-years6,  fe i(`FE') vce(cluster countyid_num)
}

foreach code in "1487" "7870" "7871" {	
	xtreg lexpend_`code' firsth prcp tmean tmax tmin months1-months6 years1-years6,  fe i(`FE') vce(cluster countyid_num)
	xtreg lexpend_`code' firsth firsth_inter_hipctblack1 prcp tmean tmax tmin months1-months6 years1-years6,  fe i(`FE') vce(cluster countyid_num)
	xtreg lexpend_`code' firsth firsth_inter_himedianhhinc prcp tmean tmax tmin months1-months6 years1-years6,  fe i(`FE') vce(cluster countyid_num)
	xtreg lexpend_`code' firsth firsth_inter_hipctbach prcp tmean tmax tmin months1-months6 years1-years6,  fe i(`FE') vce(cluster countyid_num)
	xtreg lexpend_`code' firsth firsth_inter_coastal_county prcp tmean tmax tmin months1-months6 years1-years6,  fe i(`FE') vce(cluster countyid_num)
}

restore

}

* MAIN ANALYSIS

* keep just hurricane season

keep if inrange(month,6,11)
tab month, gen(months)


if `level'==0 {

* period 2 main regressions

foreach code in "1487" "7870" "7871" {
	xtreg lexpend_`code' htreat_before hlandfall htreat_after prcp tmean tmax tmin months1-months5 years1-years6,  fe i(`FE') vce(cluster countyid_num)
}

* period2 regressions with interactions. with forecast, weather, simplified treatment definitions.

foreach code in "1487" "7870" "7871" {
	xtreg lexpend_`code' htreat_before b_inter_hipctblack1 hlandfall d_inter_hipctblack1 htreat_after a_inter_hipctblack1 prcp tmean tmax tmin months1-months5 years1-years6,  fe i(`FE') vce(cluster countyid_num)
	xtreg lexpend_`code' htreat_before b_inter_himedianhhinc hlandfall d_inter_himedianhhinc htreat_after a_inter_himedianhhinc prcp tmean tmax tmin months1-months5 years1-years6,  fe i(`FE') vce(cluster countyid_num)
	xtreg lexpend_`code' htreat_before b_inter_hipctbach hlandfall d_inter_hipctbach htreat_after a_inter_hipctbach prcp tmean tmax tmin months1-months5 years1-years6,  fe i(`FE') vce(cluster countyid_num)
	xtreg lexpend_`code' htreat_before b_inter_coastal_county hlandfall d_inter_coastal_county htreat_after a_inter_coastal_county prcp tmean tmax tmin months1-months5 years1-years6,  fe i(`FE') vce(cluster countyid_num)
	xtreg lexpend_`code' htreat_before b_inter_fema_county_week hlandfall d_inter_fema_county_week htreat_after a_inter_fema_county_week fema_county_week prcp tmean tmax tmin months1-months5 years1-years6,  fe i(`FE') vce(cluster countyid_num)
}


* period 2 analysis by common sizes. note flashlights are dominated by single counts, so doesn't make sense there.

foreach stub in "w16906" "w16912" "w128001" "b2ct" "b4ct" "b8ct" {
	xtreg lexpend_`stub' htreat_before hlandfall htreat_after prcp tmean tmax tmin months1-months5 years1-years6,  fe i(`FE') vce(cluster countyid_num)
}

* period 2 analysis by evacuation status

foreach code in "1487" "7870" "7871" {
foreach status in "evac" "nonevac" {
	xtreg lexpend_`code' htreat_before_`status' hlandfall_`status' htreat_after_`status' prcp tmean tmax tmin months1-months5 years1-years6,  fe i(`FE') vce(cluster countyid_num)
	xtreg lexpend_`code' htreat_before_`status' b_inter_coastal_county_`status' hlandfall_`status' d_inter_coastal_county_`status' htreat_after_`status' a_inter_coastal_county_`status' prcp tmean tmax tmin months1-months5 years1-years6,  fe i(`FE') vce(cluster countyid_num)
}

}

* threatened near miss vs. threatened followed by strike.

	* is this county week impacted by hurricane this week, next week, or week after that?
	egen temp2=rowmax(hlandfall hlandfall_f1 hlandfall_f2)
	* is this county week threatened by a hurricane this week, but not impacted by hurricane this week, next week, or week after that?
	gen nearmiss = (temp2==0 & htreat_before==1)
	* is this county week threatened by a hurricane this week, and subsequently impacted by hurricane this week, next week, or week after that?
	gen threatened_hit = (temp2==1 & htreat_before==1)

	* how many county weeks are hit by a storm but not threatened?
	bys countyid_num year: egen temp3=max(threatened_hit)
	count if hlandfall==1 & temp3==0
	
	* summary stats
	tab htreat_before
	tab nearmiss
	tab threatened_hit

	* generate impact leads
	tsset countyid_num week 
	gen nearmiss_l1=l1.nearmiss
	gen nearmiss_l2=l2.nearmiss

	bys week: egen temp33=max(hlandfall)
	bys week: egen temp44=max(htreat_after)

	* generated simulated impacts
	gen nearmiss_during = (nearmiss==1 | nearmiss_l1==1 | nearmiss_l2==1) & (temp33==1)
	gen nearmiss_after = (nearmiss==1 | nearmiss_l1==1 | nearmiss_l2==1) & (temp44==1)

	foreach code in "1487" "7870" "7871" {

		xtreg lexpend_`code' htreat_before hlandfall htreat_after prcp tmean tmax tmin months1-months5 years1-years6,  fe i(`FE') vce(cluster countyid_num)
		xtreg lexpend_`code' nearmiss nearmiss_during nearmiss_after prcp tmean tmax tmin months1-months5 years1-years6 ,  fe i(`FE') vce(cluster countyid_num)
		xtreg lexpend_`code' threatened_hit hlandfall htreat_after prcp tmean tmax tmin months1-months5 years1-years6 ,  fe i(`FE') vce(cluster countyid_num)
	}

}

else {

*store-level main regression
foreach code in "1487" "7870" "7871" {
	xtreg lexpend_`code' htreat_before hlandfall htreat_after prcp tmean tmax tmin months1-months5 years1-years6,  fe i(`FE') vce(cluster countyid_num)
}
}

STOP
STOP
STOP
STOP



